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Chapter 7 

Numerical methods for computing 

Casimir interactions 



Steven G. Johnson 



Abstract 

We review several different approaches for computing Casimir forces and 
^"^ related fluctuation-induced interactions between bodies of arbitrary shapes 

^^ and materials. The relationships between this problem and well known com- 

r"| putational techniques from classical electromagnet ism are emphasized. We 

Qh also review the basic principles of standard computational methods, catego- 

H_^ rizing them according to three criteria — choice of problem, basis, and solution 

5 technique — that can be used to classify proposals for the Casimir problem as 

well. In this way, mature classical methods can be exploited to model Casimir 



qH physics, with a few important modifications. 



7.1 Introduction 



Thanks to the ubiquity of powerful, general-purpose computers, large-scale 
numerical calculations have become an important part of every field of sci- 
r^ ence and engineering, enabling quantitative predictions, analysis, and de- 

sign of ever more complex systems. There are a wide variety of different 
approaches to such calculations, and there is no single "best" method for all 
•• circumstances — not only are some methods better suited to particular situa- 

. J^ tions than to others, but there are also often severe trade-offs between gen- 

^^ erality/simplicity and theoretical efficiency. Even in relatively mature areas 

^ like computational classical electromagnetism (EM), a variety of techniques 

spanning a broad range of sophistication and generality remain in widespread 
use (and new variations are continually developed) (l}|8]. Semi-analytical ap- 
proaches also remain important, especially perturbative techniques to decom- 
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pose problems containing widely differing length scales (the most challenging 
situation for brute- force numerics). Nevertheless, many commonalities and 
guiding principles can be identified that seem to apply to a range of numeri- 
cal techniques. 

Until a few years ago, Casimir forces and other EM fiuctuation-induced in- 
teractions occupied an unusual position in this tableau. Realistic, general nu- 
merical methods to solve for Casimir forces were simply unavailable; solutions 
were limited to special high-symmetry geometries (and often to special mate- 
rials like perfect metals) that are amenable to analytical and semi-analytical 
approaches. This is not to say that there were not, in principle^ decades- 
old theoretical frameworks capable of describing fiuctuations for arbitrary 
geometries and materials, but practical techniques for evaluating these the- 
oretical descriptions on a computer have only been demonstrated in the last 
few years [9]j27]. In almost all cases, these approaches turn out to be closely 
related to computational methods from classical EM, which is fortunate be- 
cause it means that Casimir computations can exploit decades of progress 
in computational classical EM once the relationship between the problems 
becomes clear. The long delay in developing numerical methods for Casimir 
interactions, from the time the phenomenon was first proposed in 1948 | 28] , 
can be explained by three factors. First, accurate measurements of Casimir 
forces were first reported only in 1997 |29 and experimental interest in com- 
plex Casimir geometries and materials has only recently experienced dramatic 
growth due to the progress in fabricating nanoscale mechanical devices. Sec- 
ond, even the simplest numerical prediction of a single force requires the 
equivalent of a large number of classical EM simulations, a barrier to casual 
numerical experimentation. Third, there have historically been many equiv- 
alent theoretical formulations of Casimir forces, but some formulations are 
much more amenable to computational solution than others, and these for- 
mulations are often couched in a language that is opaque to researchers from 
classical computational EM. 

This purpose of this review is to survey the available and proposed numer- 
ical techniques for evaluating Casimir forces, energies, torques, and related 
interactions, emphasizing their relationships to standard classical-EM meth- 
ods. Our goal is not to identify a "best" method, but rather to illuminate the 
strengths and weaknesses of each approach, highlighting the conclusions that 
can be gleaned from the classical experience. We will review an intellectual 
framework in which to evaluate different numerical techniques, comparing 
them along several axes for which quasi-independent choices of approach can 
be made. We will also emphasize a few key departures of Casimir problems 
from ordinary classical EM, such as the necessity of imaginary- or complex- 
frequency solutions of Maxwell's equations and the need for wide-bandwidth 
analyses, that impact the adaptation of off-the-shelf computational methods. 
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7.2 Characterization of numerical methods: Three axes 

Numerical methods from distinct groups or research papers often differ in 
several ways simultaneously, complicating the task of directly comparing or 
even describing them. In order to organize one's understanding of numerical 
approaches, it is useful to break them down along three axes of comparison, 
representing (roughly) independent choices in the design of a method: 

• What problem does the method solve — even within a single area such as 
classical EM, there are several conceptually different questions that one 
can ask and several ways of asking them that lead to different categories 
of methods. 

• What basis is used to express the unknowns — how the infinite number 
of unknowns in the exact partial differential equation (PDE) or integral 
equation are reduced to a finite number of unknowns for solution on a 
computer. 

• What solution technique is used to determine these unknowns — even 
with the same equations and the same unknowns, there are vast differ- 
ences among the types of direct, sparse, and iterative methods that can 
be used to attack the problem, and the efficient application of a partic- 
ular solution technique to a particular problem is sometimes a research 
task unto itself. 

In this section, we briefly summarize the available problems, basis choices, 
and solution techniques for Casimir problems. In subsequent sections, we 
then discuss in more detail the specific approaches that have currently been 
demonstrated or proposed. 



7.2.1 Posing Casimir problems 

In classical EM, there are several types of problems that are typically posed IgI 
appendix D], such as computing source-free time-harmonic eigensolutions 
E, H ~ g-*u;^ Q^j^^j eigenfrequencies cj, computing time-harmonic fields re- 
sulting from a time-harmonic current source J ~ e"*'^^, or computing the 
time-dependent fields created by an arbitrary time-dependent source J(t) 
starting at t = 0. Although these are all closely mathematically related, and 
in some sense the solution of one problem can give solutions to the other 
problems, they lead to very different types of numerical simulations. 

In a similar way, despite the fact that different formulations of Casimir- 
interaction problems are ultimately mathematically equivalent (although the 
equivalencies are often far from obvious) — and are usually answering the same 
conceptual question, such as what is the force or interaction energy for some 
geometry — each one leads most naturally to distinct classes of computational 
methods. Here, we exclude formulations such as proximity- force ("parallel- 
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plate") approximations [3Q}|32] , pairwise summation of Casimir-Polder forces 
(valid in the dilute-gas limit) (33||35) , and ray optics |36||39] , that are useful 
in special cases but represent uncontrolled approximations if they are applied 
to arbitrary geometries. Although at some point the distinctions are blurred 
by the mathematical equivalencies, we can crudely categorize the approaches 
as: 

• Computing the eigenfreq uenc ies uJn and summing the zero-point energy 

^^^ 1 28, 40. See Sec. 



7.3 



Integrating the mean energy density or force density (stress tensor), by 
evaluating field correlation functions (EiEj)^ and (HiHj)^ in terms of 
the classical EM Green's functions at cu via the fluctuation-dissipation 



theorem [Tl}|l3)|l§[20||23|[2^ See Sec. [73j 

Evaluating a path-integral expression for the interaction energy (or its 
derivative), constrained by the boundary conditions — usually, portions 
of the path integrals are performed analytically to express the problem 
in terms of classical scattering matrices or Green's functions at each uj [9J 



10,14 18,21,22,24 26 . See Sec. 7.6 



In each case, the result must be summed/integrated over all frequencies uj to 
obtain the physical result (corresponding to thermodynamic/quantum fluc- 
tuations at all frequencies). The relationship of the problem to causal Green's 
functions (fields appear after currents) means that the integrand is analytic 
for Imcj > |41 . As a consequence, there is a choice of contours of u inte- 
gration in the upper-half complex plane, which is surprisingly important — it 
turns out that the integrands are wildly oscillatory on the real-cj axis and re- 
quire accurate integration over a huge bandwidth, whereas the integrands are 
much better-behaved along the imaginary-cj axis ("Wick-rotated" or "Matsub- 
ara" frequencies) . This means that Casimir calculations almost always involve 
classical EM problems evaluated at complex or imaginary frequencies^ as is 
discussed further in Sec. |7.4| The nonzero-temperature case, where the inte- 
gral over imaginary frequencies becomes a sum (numerically equivalent to a 
trapezoidal-rule approximation), is discussed in Sec. |7.8[ 

There is also another way to categorize the problem to be solved: whether 
one is solving a partial differential equation (PDE) or an integral equa- 
tion. In a PDE, one has volumetric unknowns: fields or other functions at 
every point in space, related to one another locally by derivatives and so on. 
In an integral equation, one typically has surface unknowns: the fields or 
currents on the boundaries between piecewise- homogeneous regions, related 
to one another non-locally by the Green's functions of the homogeneous re- 



gions (typically known analytically) [ll^ (described further in Sec. 7.5.3). 
The key point is to take advantage of the common situation in which one has 
piecewise-constant materials, yielding a surface integral equation. (There are 
also volume integral equations for inhomogeneous media [42j, as well as hybrid 
integral/PDE approaches |1 , but these are less common.) There are other 



hybrid approaches such as eigenmode expansion 143^-45 , also called rigorous 
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coupled-wave analysis (RCWA) | 46p7l or a cross-section method (48) : a struc- 
ture is broken up along one direction into piecewise-constant cross-sections, 
and the unknown fields at the interfaces between cross-sections are propa- 
gated in the uniform sections via the eigenmodes of those cross-sections (com- 
puted analytically or numerically by solving the PDE in the cross-section). 
Eigenmode expansion is most advantageous for geometries in which the cross- 
section is constant over substantial regions, just as integral-equation methods 
are most advantageous to exploit large homogeneous regions. 



7.2.2 Choices of basis 

Casimir problems, for the most part, reduce to solving classical EM lin- 
ear PDEs or integral equations where the unknowns reside in an infinite- 
dimensional vector space of functions. To discretize the problem approxi- 
mately into a finite number N of unknowns, these unknown functions must 
be expanded in some finite basis (that converges to the exact solution as 
N -^ oo). There are three typical types of basis: 



• Finite differences Pl|^^^ (FD): approximate a function f{x) by 
its values on some uniform grid with spacing Z\x, approximate deriva- 
tives by some difference expression [e.g. second-order center differences 
f'{x) ~ ^ — ^2Ax ^ -\-0{Ax'^)] and integrals by summations (e.g. a 
trapezoidal rule). 

• Finite-element methods jT^'s','?,'?! (FEM): divide space into geometric 
elements (e.g. triangles/tetrahedra), and expand an unknown f{x) in a 
simple localized basis expansion for each element (typically, low-degree 
polynomials) with some continuity constraints. (FD methods are viewable 
as special cases of FEMs for uniform grids.) For an integral-equation 
approach, where the unknowns are functions on surfaces, the same idea 
is typically called a boundary-element method (BEM) [Tpl[7l5l|52| p] 

• Spectral methods |53 : expand functions in a non- localized complete 
basis, truncated to a finite number of terms. Most commonly, Fourier 
series or related expansions are used (cosine series, Fourier-Bessel series, 
spherical or spheroidal harmonics, Chebyshev polynomials, etc.). 

Finite differences have the advantage of simplicity of implementation and 
analysis, and the disadvantages of uniform spatial resolution and relatively 
low-order convergence (errors typically ~ Ax'^ |2 or even ~ Ax in the pres- 
ence of discontinuous materials unless special techniques are used (54|[55]). 
FEMs can have nonuniform spatial resolution to resolve disparate feature 



^ The name method of moments is also commonly applied to BEM techniques for EM. How- 
ever, this terminology is somewhat ambiguous, and can refer more generally to Galerkin 
or other weighted-residual methods (and historically referred to monomial test functions, 
yielding statistical "moments") 53j. 



6 Steven G. Johnson 

sizes in the same problem, at a price of much greater complexity of imple- 
mentation and solution techniques, and can have high-order convergence at 
the price of using complicated curved elements and high-order basis func- 
tions. Spectral methods can have very high-order or possibly exponential 
("spectral") convergence rates ^3] that can even suit them to analytical 
solution — hence, spectral methods were the dominant technique before the 
computer era and are typically the first class of methods that appear in any 
field, such as in Mie's classic solution of wave scattering from a sphere |56 . 
However, exponential convergence is usually obtained only if all disconti- 
nuities and singularities are taken explicitly into account in the basis |53 . 
With discontinuous materials, this is typically only practical for very smooth, 
high-symmetry geometries like spheres, cylinders, and so on; the use of a 
generic Fourier/spectral basis for arbitrary geometries reduces to a brute- 
force method that is sometimes very convenient |57|, but may have unremark- 
able convergence rates [53l[57{[58] . BEMs require the most complicated imple- 
mentation techniques, because any nontrivial change to the Green's functions 
of the homogeneous regions (e.g. a change in dimensionality, boundary con- 
ditions, or material types) involves tricky changes to the singular-integration 
methods required to assemble the matrix |59-61 and to the fast-solver meth- 
ods mentioned in Sec. 17.2.31 

Given FEM/BEM or spectral basis functions bn{x) and a linear equa- 
tion Au{x) = v{x) for an unknown function u in terms of a linear differen- 
tial/integral operator A, there are two common ways |53 to obtain a finite 
set of N equations to determine the N unknown coefficients c^ in u{x) ~ 
Xln^^^^(^)- ^^^ ^^ ^ collocation method: require that {Au — v)\x^ = be 
satisfied at N collocation points Xn- The other is a Galerkin method: require 
that (6/e, Au — 'u) = be satisfied for /c = 1, . . . , A^, where (•, •) is some inner 
product on the function space. Both approaches result in an A^ x A^ matrix 
equation of the form Au = v. A Galerkin method has the useful property 
that if A is Hermitian and/or definite then the matrix A^n = {^ki^^n) has 
the same properties. 

The specific situation of vector- valued unknowns in EM creates additional 
considerations for the basis functions. In order to obtain center-difference 
approximations for all the field components, FD methods for EM typically 
use a staggered Yee grid 1 2 , 49 , in which each component of the EM fields 
is offset onto its own ^-shifted grid. In FEMs for EM, in order to main- 
tain the appropriate continuity conditions for curl or divergence operators, 
one uses special classes of vector-valued basis functions such as Ni/^oedi'/^oelec 
elements (7||62]. In BEMs for EM, vector- valued RWG (Rao, Wilton, and 
Glisson) basis functions [63] (or generalizations thereof |64 ) are used in or- 
der to enforce a physical continuity condition on surface currents (to preclude 
accumulation of charge at element edges); see also Fig. 7.3 in Sec. |7.5.3 A 



spectral integral-equation method for EM with cylindrical or spherical scat- 
terers is sometimes called a multipole-expansion method |5|, since the 
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obvious spectral basis is equivalent to expanding the scattered fields in terms 
of multipole moments. 



7.2.3 Solution techniques for linear equations 

Given a particular problem and basis choice, one at the end obtains some N x 
N set of linear equations Ax = b to solve (or possibly eigenequations Ax = 
A5x)j^Note also that a single Casimir-force calculation requires the solution 
of many such equations, at the very least for an integral over frequencies (see 
Sec. |7.4[). There are essentially three ways to solve such a set of equations: 



• Dense-direct solvers: solve Ax = b using direct matrix-factorization 
methods (e.g. Gaussian elimination) rlrequiring 0{N'^) storage and 0{N^) 
time |65 . 

• Sparse-direct solvers i66|: if A is sparse (mostly zero entries), use sim- 
ilar direct matrix- factorization methods, but cleverly re- arranged in an 
attempt to preserve the sparsity. Time and storage depend strongly on 
the sparsity pattern of A (the pattern of nonzero entries). 

• Iterative methods l65J[67Jl68 : repeatedly improve a guess for the solu- 
tion X (usually starting with a random or zero guess), only referencing 
A via repeated matrix-vector multiplies. Time depends strongly on the 
properties of A and the iterative technique, but typically requires only 
0{N) storage. Exploits any fast way [ideally 0{N) or 0{N log N)] to 
multiply A by any arbitrary vector. 

If the number N of degrees of freedom is small, i.e. if the basis converges 
rapidly for a given geometry, dense-direct methods are simple, quick, and 
headache-free (and have a standard state-of-the-art implementation in the 
free LAPACK library [69]). For example, N = 1000 problems can be solved 
in under a second on any modern computer with a few megabytes of mem- 
ory. Up to A^ ~ 10^ is reasonably feasible, but N = 10^ requires almost 
100 GB of memory and days of computation time without a large parallel 
computer. This makes dense-direct solvers the method of choice in simple 
geometries with a rapidly converging spectral basis, or with BEM integral- 
equation methods for basic shapes that can be accurately described by a few 
thousand triangular panels, but they rapidly become impractical for larger 
problems involving many and/or complex objects (or for moderate-size PDE 
problems even in two dimensions). 



^ This applies equally well, if somewhat indirectly, to the path-integral expressions of 
Sec. |7.6| where one evaluates a log determinant or a trace of an inverse, since this is done 
using either eigenvalues or the same matrix factorizations that are used to solve Ax = b. 
^ Technically, all eigensolvers for A^ > 4 are necessarily iterative, but modern dense- 
eigensolver techniques employ direct factorizations as steps of the process [65]. 
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In PDE methods with a locahzed (FD or FEM) basis, the matrices A 
have a special property: they are sparse (mostly zero). The locality of the 
operators in a typical PDE means that each grid point or element directly 
interacts only with a bounded number of neighbors, in which case A has only 
0{N) nonzero entries and can be stored with 0{N) memory. The process of 
solving Ax = b, e.g. computing the LU factorization A = LU by Gaussian 
elimination |65 , unfortunately, ordinarily destroys this sparsity: the resulting 
L and U triangular matrices are generally not sparse. However, the pattern 
of nonzero entries that arises from a PDE is not random, and it turns out 
that clever re-orderings of the rows and columns during factorization can par- 
tially preserve sparsity for typical patterns; this insight leads to sparse- direct 
solvers |66 , available via many free-software packages implementing different 
sparsity-preserving heuristics and other variations |68 . The sparsity pattern 
of A depends on the dimensionality of the problem, which determines the 
number of neighbors a given element interacts with. For meshes/grids having 
nearest-neighbor interactions, a sparse-direct solver typically requires 0{N) 
time and storage in Id (where the matrices are band-diagaonal), 0{N^^'^) 
time with 0{N log N) storage in 2d, and 0{N'^) time with 0(A^^/^) storage 
in 3d (GGJITO]. The practical upshot is that sparse-direct methods work well 
for Id and 2d PDEs, but can grow to be impractical in 3d. For BEM and 
spectral methods, the interactions are not localized and the matrices are not 
sparse, so sparse-direct methods are not directly applicable (but see below 
for an indirect technique). 

For the largest-scale problems, or for problems lacking a sparse A, the re- 
maining possibility is an iterative method. In these methods, one need only 
supply a fast way to multiply A by an arbitrary vector y, and the trick is to 
use this Ay operation on a clever sequence of vectors in such a way as to make 
an arbitrary initial guess xq converge as rapidly as possible to the solution 
X, ideally using only 0{N) storage. Many such techniques have been devel- 



oped 165,67,68 . The most favorable situation for Ax = b occurs when A is 
Hermit ian positive-definite, in which case an ideal Krylov method called the 
conjugate- gradient method can be applied, with excellent guaranteed con- 
vergence properties [65n67j, and fortunately this is precisely the case that 
usually arises for the imaginary-frequency Casimir methods below. There are 
two wrinkles that require special attention, however. First, one must have a 
fast way to compute Ay. If A is sparse (as for PDE and FD methods), then 
only 0{N) nonzero entries of A need be stored (as above) and Ay can be 
computed in 0{N) operations. In a spectral method, A is generally dense, 
but for spectral PDE methods there are often fast 0{NlogN) techniques to 
compute Ay using only 0{N) storage {A is stored implicitly), based on fast 
Fourier transform (FFT) algorithms (53J|57]. In a BEM, where A is again 
dense, a variety of sophisticated methods that require only 0(A^log A^) com- 
putation time and 0{N) storage to compute Ay (again storing A implicitly) 
have been developed |1,3,7,71 , beginning with the pioneering fast-multipole 
method (FMM) 172J. These fast BEMs exploit the localized basis and the 
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decaying, convolutional nature of the Green's function to approximate long- 
range interactions (to any desired accuracy). FMMs can be viewed as an 
approximate factorizations into sparse matrices, at which point sparse-direct 
methods are also applicable |73|. A second wrinkle is that the convergence 
rates of iterative methods depend on the condition number of A (the ratio of 
largest and smallest singular values) |65[|67|, and condition numbers generally 
worsen as the ratio of the largest and smallest lengthscales in the problem 
increases. To combat this, users of iterative methods employ precondition- 
ing techniques: instead of solving Ax = 6, one solves KAx = Kb or similar, 
where the preconditioning matrix K is some crude approximate inverse for 
A (but much simpler to compute than A~^\) such that the condition num- 
ber of KA is reduced |67 . The difficulty with this approach is that good 
preconditioners tend to be highly problem-dependent, although a variety of 
useful approaches such as incomplete factorization and coarse-grid/multigrid 
approximations have been identified i65ii67i. The upshot is that, while the 
largest-scale solvers almost invariably use iterative techniques, for any given 
class of physical problems it sometimes takes significant research before the 
iterative approach becomes well-optimized. 



7.3 The impracticality of eigenmode summations 

Perhaps the simplest way to express the Casimir energy, at zero temperature, 
is as a sum of zero-point energies of all oscillating EM modes in the system: 



U'E"^. P.i) 



where uj is the frequency of the mode (~ e"*'^^) [28l[74|. That is, when the 
electromagnetic field is quantized into photons with energy hw^ it turns out 
that the vacuum state in the absence of photons is not empty, but rather 
has the energy equivalent of "half a photon" in each mode. The computa- 
tional strategy is then straightforward, in principle: compute the EM eigen- 
frequencies uj in the problem by some numerical method (many techniques 
are available for computing eigenfrequencies |1,57 ) and sum them to obtain 
U. Forces are then given by the derivative of U with respect to changes in the 
geometry, which could be approximated by finite differences or differentiated 
analytically with a Hellman-Feynman technique 1 7^ (more generally, deriva- 
tives of any computed quantity can be computed efficiently by an adjoint 
method 176 ). 



Of course, U in eq. (7.1) has the disadvantage of being formally infinite, 
but this is actually a minor problem in practice: as soon as one discretizes the 
problem into a finite number of degrees of freedom (e.g., a finite number of 
grid points) , the number of eigenfrequencies becomes finite (with the upper 
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bound representing a Nyquist-like frequency of the grid). This is the nu- 
merical analogue |12 of analytical regularization techniques that are applied 
to truncate the same sum in analytical computations f28^. (These regular- 
izations do not affect energy differences or forces for rigid-body motions.) 
Matters are also somewhat subtle for dissipative or open systems |77|. But 
the most serious problem is that, even in the lossless case, this sum is badly 
behaved: even when one differentiates with separation a to obtain a finite 
force F = — f XI ^ ? ^^^ summand is wildly oscillatory and includes sub- 
stantial contributions from essentially every frequency, which mostly cancel 
to leave a tiny result |T2J|78]. Numerically, therefore, one must ostensibly 
compute all of the modes, to high precision, which requires 0{N^) time and 
0{N'^) storage (for a dense-direct eigensolve r \65{) given N degrees of free- 
dom. This is possible in simple Id problems [12,40 , but is impractical as a 
general approach. 

Because of the mathematical equivalence of the different approaches to 
the Casimir problem, the mode-summation method is sometimes useful as 
a starting point to derive alternative formulations, but the end result is in- 
variably quite different in spirit from computing the eigenfrequencies one by 
one and summing them. For example, if one has a function z{uj) whose roots 
are the eigenfrequencies, then one can equivalently write [/, via the residue 
theorem of complex analysis, as /7 = ^ ^^ %^ doj^ ^'^^ where C is any 
closed contour in the complex-cj plane that encloses the roots [79 1. However, 
finding functions whose roots are the eigenfrequencies naturally points to- 
wards Green's functions (to relate different boundary conditions), and the 
contour choices typically involve Wick rotation as in Sec. |7.4[ so this ap- 
proach leads directly to imaginary-frequency scattering-matrix techniques as 
in Sec. 7.6 |18 . A similar contour integral arises from a zeta-function regu- 



larization of ([71]) [80j. 



7.4 The complex-frequency plane and contour choices 

In order to better understand the frequency integration/summation in Casimir 
problems, it is illustrative to examine the analytical formula for the simple 
case of two perfect-metal plates in vacuum separated by a distance a, in which 
case it can be derived in a variety of ways that the attractive force F is given 
by (81]: 
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h 


3 Re 


/ dco 

Jo 


/"CO 

/ * 


Re 


" />oo 

7o 


f{uj)duj 


= Im 



p^uj^ 



o2ip{uj+i0+)a/c 



fiiOd^ 



he 

240a4 ' 



(7.2) 



7 Numerical methods for computing Casimir interactions 11 

where f{uj) is the contribution of each frequency u to the force and p is re- 
lated to the plate-parallel momentum of the contributing modes/fluctuations. 
In this special case, the entire integral can be performed analytically, but for 
parallel plates of some finite permittivity £ the generalization (the Lifshitz 
formula |81 ) must be integrated numerically. In practice, however, the for- 
mula and its generalizations are never integrated in the form at left — instead, 
one uses the technique of contour integration from complex analysis to Wick 
rotate the integral to imaginary frequencies uo = i^^ integrating over ^. (In 
fact, the formula is typically derived starting in imaginary frequencies, via 
a Matsubara approach |81 .) In this section, we review why a trick of this 
sort is both possible and essential in numerical computations for all of the 
methods described below. 

Wick rotation is always possible as a consequence of causality. It turns out 
that the frequency contributions /(cj) for arbitrary materials and geometries, 
for all of the different formulations of the Casimir force below, are ultimately 
expressed in terms of classical EM Green's functions at uo: the EM fields 
in response to time-harmonic currents J r^ e~^^'^ . As a consequence of the 
causality of Maxwell's equations and physical materials — EM fields always 
arise after the source currents, not before — it mathematically follows that 
the Green's functions must be analytic functions (no poles or other singu- 
larities) when Imcj > (the upper-half complex plane) |4l|. Poles in the 
Green's function correspond to eigenfrequencies or resonances of the source- 
free Maxwell's equations, and must lie at Imo; < for any physical system 
with dissipative materials (with the poles approaching Imcj = 0~ in the ide- 
alized lossless limit). [One can easily see explicitly that this is true for the 
f{uo) above: the poles result from a vanishing denominator in the p integrand, 
which only occurs for purely real uj corresponding to the real-frequency modes 
trapped between two perfect-metal plates.] As an elementary consequence of 
complex analysis, this analyticity means that the J duo can be arbitrarily de- 
formed to any contour in the upper-half complex-cj plane without changing 
the integration result. 

Wick rotation is essential for computation because the frequency contri- 
butions f{(jj) to the force (or interaction energy or other related quantities) 
are extremely ill-behaved near to the real-cj axis: they are wildly oscillatory 
and slowly decaying. For example, the magnitude and phase of the function 



f{uj) are plotted in the complex uj plane in Fig. 7.1, where the p integral 
was evaluated numerically with a high-order Clenshaw-Curtis quadrature 
scheme |82 . Merely evaluating f{uo) along the real-cj axis is difficult because 
of singularities (which ultimately reduce the integral to a summation over 



eigenfrequency-contributions as in Sec. 7.3); in physical materials with dissi- 
pation, the real-cj axis is non-singular but is still badly behaved because of 
poles (lossy modes) located just below the axis. Along any contour parallel 
to the real-cj axis, the integrand is oscillatory (as can be seen from the phase 
plot) and non-decaying (as can be seen from the magnitude plot): formally, 
just as with the infinite summation over eigenmodes in Sec. |7.3[ one must 
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Fig. 7.1 Contributions f{uj) to the Casimir force, from each fluctuation/mode frequency 
a;, for two perfect-metal plates with separation a, in the complex-o; plane. Left: magnitude 
\f{uj)\. Right: phase Zf(uj). [The magnitude is truncated at 10^/i/a^, as it diverges towards 
the real-a; axis, and some numerical artifacts (rapid oscillations) are visible near the real-o; 
axis in the phase due to difficulty in evaluating fioj). ] The key point is that f{(jj) is badly 
behaved (oscillatory and non-decaying) along contours parallel to the real-o; axis, whereas 
/(a;) is nicely behaved (non-oscillatory and exponentially decaying) along contours parallel 
to the imaginary-o; axis. 



integrate over an infinite bandwidth, regularized in some way (e.g. by the 
Nyquist frequency placing an upper bound on uo for a finite grid), where the 
oscillations almost entirely cancel to leave a tiny remainder (the force). (Any 
physical materials must cease to polarize as cj ^ oo where the susceptibility 
vanishes |41 , which will make the force contributions eventually vanish as 
a; ^ oo even in Id, but a very wide-bandwidth oscillatory integral is still 
required.) This is a disaster for any numerical method — even when one is 
only integrating an analytical expression such as the Lifshitz formula, mere 
roundoff errors are a severe difficulty for real uj. Along the imaginary-cj axis, 
on the other hand (or any sufficiently vertical contour), f{uo) is exponen- 
tially decaying and mostly non-oscillatory — an ideal situation for numerical 
integration. 

Therefore, in order for classical EM solvers to be used for Casimir prob- 
lems, they need to be adapted to solve Maxwell's equations at complex or 
imaginary uj. Although this sounds strange at first, the frequency-domain 
problem actually becomes numerically easier in every way at imaginary a;; 
this is discussed in more detail in Sec. |7.5.1.2| In fact, one can even identify 
an exact mathematical equivalence between a particular complex-cj contour 
and a rea/-frequency system where an artificial dissipation has been intro- 
duced, as discussed in Sec. |7.5.5| below — using this trick, one can actually 
use classical EM solvers with no modification at all, as long as they han- 
dle dissipative media. In any case, one needs an integral over frequencies to 
compute a physically meaningful quantity, which means that solvers and ma- 
terial models, not to mention any physical intuition used for guidance, must 
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be valid for more than just a narrow real-o; bandwidth (unlike most problems 
in classical EM). 

Numerically, it should be pointed out that the f{i£,) integrand is smooth 
and exponentially decaying, and so the ^ integral can be approximated to high 
accuracy by an exponentially convergent quadrature (numerical integration) 
scheme using evaluations at relatively few points ^. For example, one can use 
Gauss-Laguerre quadrature |83 , Gaussian quadrature with an appropriate 
change of variables |84 , or Clenshaw-Curtis quadrature with an appropriate 
change of variables 1 82 . 



7.5 Mean energy/force densities and the 
fluctuation— dissipation theorem 

Another, equivalent, viewpoint on Casimir interactions is that they arise 
from geometry-dependent fluctuations of the electromagnetic fields E and 
H, which on average have some nonzero energy density and exert a force. If 
we can compute these average fields, we can integrate the resulting energy 
density, stress tensors, and so on, to obtain the energy, force, or other quan- 
tities of interest. The good news is that there is a simple expression for those 
fluctuations in terms of the fluctuation— dissipation theorem of statistical 
physics: the correlation function of the fields is related to the corresponding 
classical Green's function |81 . Ultimately, this means that any standard clas- 
sical EM technique to compute Green's functions (fields from currents) can 
be applied to compute Casimir forces, with the caveat that the techniques 
must be slightly modified to work at imaginary or complex frequencies as 
described below. 



7.5.1 Background 

The temperature-T correlation function for the fluctuating electric field at a 
given frequency u is given by [sT] : 

{Ej{^)Ek{^'))^ = --Im [uj^Gff^iuj;^,^')] coth{huj/2k^T), (7.3) 

where G^ = {G^)j is the classical dyadic "photon" Green's function, pro- 
portiona]^ to the relationship between an electric-dipole current in the k 
direction at x' to the electric field at x, and solves 



"^ The electric field E(x) from a dipole current J = (5^(x — x^je^e ^'^^ is E(x) 
ia;Gf(a;,x,xOe-^^*. 
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[V X /i(a;, x)"^ V x - uj'^e{uj, x)] Gf (cj, x, x') = (5^(x - x')e/c, (7.4) 

where e is the electric permittivity tensor, fi is the magnetic permeabihty 
tensor, and e^ is a unit vector in direction k. Similarly, the magnetic-field 
correlation function is 

{Hj{^)Hk{^'))^ = --Im [cj2Gfjcj;x,xO] coth{hw/ 2 k^T). (7.5) 

The magnetic Green's function G^ can be defined in two essentially equiva- 
lent ways. The first is as derivatives 2'^f \ ^ x G^ x V^ —t-tt of the electric 
Green's function G£,(x, x'), where V and V denote derivatives with respect 
to X and x' (V' acting to the left), respectively |81 . The second way to de- 
fine G^ is proportional to the magnetic field in response to a magnetic-dipole 



current, analogous to (7.4): 



[V X s{uj, x)-^ V X - cjV(^, x)] Gf (cj, X, xO = (^^(x - x^e/c, (7.6) 

which can be more convenient for numerical calculation ll3i. These two defi- 
nitions are related 186 by G^ = 2^ ^ V x G^ x V'^^ 2-^r^(^(x-xO/ 

(with / being the 3x3 identity matrix) jj where the second (diagonal) term 
has no effect on energy differences or forces below and is therefore irrelevant. 
Now, these equations are rather nasty along the real-o; axis: not only will 
there be poles in G just below the axis corresponding to lossy modes, but 
in the limit where the dissipative losses vanish {e and /i become real), the 
combination of the poles approaching the real axis with the Im in the corre- 
lation function results in a delta function at each polqj and integrals of the 
correlation functions turn into sums over modes as in Sec. |7.3[ However, the 
saving grace, as pointed out in Sec. |7.4[ is that Green's functions are causal^ 
allowing us to transform any integral over all real fluctuation frequencies into 
an integral over imaginary fluctuation frequencies uj = i^. The coth factor 
has poles that alter this picture, but we will eliminate those for now by con- 
sidering only the T = 0+ case where coth(+oo) = 1, returning to nonzero 



temperatures in Sec. 7.8 



^ This can be seen more exphcity by substituting G^^ = ^ly x G^^ x V^^ ^S 



into ( |7.6| ), with S denoting (5(x — x^)/ and fi or fi' denoting /x(x) or /i(x^), respectively. In 
particular, [V x ^ V x -uj'^Li\(^^Vx G^ x V A - -^S) yields V x [4- V x ^ V x G^ - 
G^l X V^^-Vx J- , \/xd + 6, which via (TmI gives +V x ^(5x V^^-Vx -J-^V x 

^ /J. uj^iJ,'e ' |____r ^^ /^ oj'^iJi'e 

5 -\- 5 = 5 diS desired, where in the last step we have used the fact that 5 x S/' = V x 5 
[since V x is antisymmetric under transposition and V^5(x — x^) = — V(5(x — x^)]. 
^ This follows from the standard identity that the limit Im[(a:: + 20+)"-*^], viewed as a 
distribution, yields —t:5{x) 85 . 
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7.5.1.1 Energy density 

In particular, to compute the Casimir energy /7, we merely integrate the clas- 
sical energy density in the EM field |41 over all positions and all fiuctuation 
frequencies, Wick-rotated to an integral over imaginary frequencies, resulting 
in the expression: 



U 



-!>l\ 



^(|EPk + ^(|H|=) 



i^ 



d^x, (7.7) 



where we have simplified to the case of isotropic materials (scalar e and 
ji). At thermodynamic equilibrium, this expression remains valid even for 
arbitrary dissipative/dispersive media thanks to a direct equivalence with 
a path-integral expression fST^, which is not obvious from the classical view- 
point in which the energy density is usually only derived in the approximation 
of negligible absorption |41 . (Thanks to the relationship between the Green's 
function and the local density of states 1 88 , there is also a direct equivalence 
between this energy integral and eigenmode summation |12 .) In the com- 
mon case where ji has negligible frequency dependence (magnetic responses 
are usually negligible at the short wavelengths where Casimir forces are im- 
portant, so that /i ?^ /io), we can use the identit^Tlthat j£|Ep = //i|Hp for 
fields at any given frequency |6 to simplify this expression to |12 : 



U 



=r*/4^<i'=i^>-''''^- '" 



Here, the zero-temperature imaginary- frequency mean-square electric field is 
given by: 

{|E(x)|2),j = VtrG^(ie;x,x), (7.9) 

TT 

where tr denotes the trace ^ • G^ and the Im has disappeared compared to 



(7.3) because G (i^) is real and the Im cancels the i m duj ^ id^. 



Equation (7.13) may at first strike one as odd, because one is evaluating 
the Green's function (field) at x from a source at x, which is formally infinite. 
This is yet another instance of the formal infinities that appear in Casimir 
problems, similar to the infinite sum over modes in Sec. |7.3[ In practice, this 
is not a problem either analytically or numerically. Analytically, one typi- 



^ Lest the application of this field identity appear too glib, we can also obtain the 
same equality directly from the Green's functions in the correlation functions. We have 



//i(|H|2) = |tr/^2^Qi/(x,x), and from the identity after ^M we know that C^/iG^ 

—V X G X V^— r + 5. However, because V x is self-adjoint 6 , we can integrate by 

/J, 

parts to move V x from the first argument/index of G^ to the second, obtaining 
— G-^ X V' ^ X V' = ^'^e'G^ — 5 from the first term under the integral. (Here, we em- 
ploy the fact that G^^ is real-symmetric at imaginary u = i^, from Sec. |7.5.1.2[ to apply 
( |7.10| ) to the second index/argument instead of the first.) This cancels the other delta from 
^^~/IG^ and leaves ^^sG^^, giving £(|Ep) as desired. 
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cally regularizes the problem by subtracting off the vacuum Green's function 
(equivalent to only looking at the portion of the fields at x which are refiected 
off of inhomogeneities in £ or /i) |81 . Numerically, in an FD or FEM method 
with a finite grid, the Green's function is everywhere finite (the grid is its 
own regularization) |12 . In a BEM, the Green's function is explicitly written 
as a sum of the vacuum field and scattered fields, so the former can again 



be subtracted analytically |l2]. As in Sec. 7.3, these regularizations do not 
affect physically observable quantities such as forces or energy differences, 
assuming rigid-body motion. 



7.5.1.2 The remarkable imaginary- frequency Green's function 

This imaginary-frequency Green's function is actually a remarkably nice ob- 



ject. Wick-rotating eq. (7.4), it satisfies: 



[V X /i(z^, x)-i V X + ee{i^^ x)] Gf (ze, x, x^ = ^^(x - x^e^. (7.10) 

Because of causality, it turns out that e and fi are strictly real-symmetric 
and positive-definite (in the absence of gain) along the imaginary-frequency 
axis, even for dissipative/dispersive materials |41 . Furthermore, the operator 
V X /i~^V X is real-symmetric positive-semidefinite for a positive-definite 



real-symmetric /i |6 . Thus, the entire bracketed operator [• • • ] in eq. (7.10) 
is real- symmetric positive- definite for ^ > 0, which lends itself to some of the 
best numerical solution techniques (Cholesky decomposition [65], tridiagonal 
QR |65 , conjugate gradients |65ji67], and Rayleigh-quotient methods \68\). 
(This definiteness is also another way of seeing the lack of poles or oscillations 
for uj = i(_.) It follows that the integral operator whose kernel is G^, i.e. the 



inverse of the [• • • ] operator in eq. (7.10), is also real-symmetric positive- 



definite, which is equally useful for integral-equation methods. 

In vacuum, the 3d real-cj Green's function ~ g«u;|x-x \/c^^^ _ -j^/| |4]| 
is Wick-rotated to ~ e~^^^~^ l/^/|x — x'|, an exponentially decaying, non- 
oscillatory function. This is yet another way of understanding why, for uj = i^, 
there are no interference effects and hence no "modes" (poles in G), and in- 
tegrands tend to be non-oscillatory and exponentially decaying (as ^ ^ oo, 
G becomes exponentially short-ranged and does not "see" the interacting ob- 
jects, cutting off the force contributions). (It also means, unfortunately, that 
a lot of the most interesting phenomena in classical EM, which stem from 
interference effects and resonances, may have very limited consequences for 
Casimir interactions.) 

One other property we should mention is that the operator becomes 
semidefinite for ^ = 0, with a nullspace encompassing any static field dis- 
tribution (V(/) for any scalar (j)). This corresponds to the well-known singu- 
larity of Maxwell's equations at zero frequency |89,90 , where the electric 
and magnetic fields decouple [41j. Since we eventually integrate over ^, the 
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Fig. 7.2 Schematic depiction of two objects whose Casimir interaction is desired. One 
computational method involves integrating a mean stress tensor around some closed surface 
(dashed red line) surrounding an object, yielding the force on that object. 



measure-zero contribution from ^ = does not actually matter, and one 
can use a quadrature scheme that avoids evaluating ^ = 0. However, in the 
nonzero-temperature case of Sec. |7.8| one obtains a sum over discrete-^ con- 
tributions, in which case the zero- frequency term is explicitly present. In this 
case, ^ = can be interpreted if necessary as the limit ^ ^ 0+ (which can be 
obtained accurately in several ways, e.g. by Richardson extrapolation |91 , al- 
though some solvers need special care to be accurate at low frequency [89 .90] ); 
note, however, that there has been some controversy about the zero- frequency 
contribution in the unphysical limit of perfect/dissipationless metals (921. 



7.5.1.3 Stress tensor 



In practice, one often wants to know the Casimir force (or torque) on an 
object rather than the energy density. In this case, instead of integrating 
an electromagnetic energy density over the volume, one can integrate an 
electromagnetic stress tensor over a surface enclosing the object in question. 



schematically depicted in Fig. 7.2 |81 . The mean stress tensor for the Casimir 

force is [81j: 



(^j/c(x))a; =e(x,Cj) 



(^,(x)^,(x))^-^^(^,(x) 



+ /i(x,cj) 



Sjk 



(i7,(x)i7,(x))-^^(i7,(x)2) 



. (7.11) 



As above, the field correlation functions are expressed in terms of the classical 
Green's function, and the integral of the contributions over all u is Wick- 
rotated to imaginary frequencies uo = i(^: 



^^ surface 



(7.12) 
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with (zero-temperature) correlation functions 

{^,(x)Efe(x))ie = -e^Gffc(ie;x,x), (7.13) 

{ff,(x)fffe(x)),5 = ^e'Gffc(ie;x,x), (7.14) 

corresponding to the fields on the stress-integration surface in response to 
currents placed on that surface. To compute a Casimir torque around an 
origin r, one instead uses (x — r) x (T(x))i^ • dS |93 . 

The derivation of this stress tensor ( |7.11[ ) is not as straightforward as it 
might at first appear. If the stress-integration surface lies entirely in vac- 



uum e ^ sq and /j. ^ /j^q^ then one can interpret (7.11) as merely the ordi- 
nary EM stress tensor from the microscopic Maxwell equations |41 , albeit 
integrated over fiuctuations. If the stress-integration surface lies in a disper- 
sive/dissipative medium such as a fiuid, however, then the classical EM stress 
tensor is well known to be problematic [41^ and (7.11 ) may seem superficially 



incorrect. However, it turns out that these problems disappear in the context 
of thermodynamic equilibrium, where a more careful free-energy derivation of 



the Casimir force from fiuctuations indeed results in equation (7.11) |81[|94| 
which has also proved consistent with experiments |95^,^. Note also that, 
while ( |7.11| ) assumes the special case of isotropic media at x, it can still be 
used to evaluate the force on objects made of anisotropic materials, as long 
as the stress-integration surface lies in an isotropic medium (e.g. vacuum or 
most fiuids). 

This formulation is especially important for methods that use an iterative 
solver for the Green's functions as discussed below, because it only requires 
solving for the response to currents on the stress-integration surface, rather 
than currents at every point in space to integrate the energy density, greatly 
reducing the number of right-hand sides to be solved for the linear equa- 



tion (7.10) 1 12 1; additional reductions in the number of right-hand sides are 



described in Sec. 17.5.61 

7.5.2 Finite- difference frequency- domain (FDFD) 

In order to determine the Casimir energy or force, one evaluates the Green's 



function by solving eq. ( 7.10 ) and then integrates the appropriate energy/force 
density over volume/surface and over the imaginary frequency ^. The cen- 
tral numerical problem is then the determination of the Green's function by 



solving a set of linear equations corresponding to ( |7.10 ), and probably the 



simplest approach is based on a finite-difference (FD) basis: space is divided 



^ If compressibility of the fluid and the density-dependence of e is not neglected, then there 
is an additional de/dp term in (|7.11|) resulting from fluctuations in the density p i81|. 
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into a uniform grid with some resolution Z\x, derivatives are turned into dif- 
ferences, and one solves (|7.10) by some method for every desired right-hand 



side. In classical EM (typically finding the field from a given current at real 
cj), this is known as a finite-difference frequency-domain (FDFD) method, 
and has been widely used for many years [5,49 



For example, in one dimension for z-directed currents/fields with /i = 1, 



(7.10) becomes (-^ ^^A sGf^ = J (x-x')i. If we approximate Gf^(nzAx, x') 
(jr'n, then the corresponding finite- difference equation, with a standard center- 



difference approximation for d? /dx^ |50 , is 



AX^ . s ^n^n ^^ 



e^£„G„ = ^-^, (7.15) 



replacing the 5{x — x') with a discrete equivalent at n' . Equation (7.15) is 
a tridiagonal system of equations for the unknowns Gn- More generally, of 
course, one has derivatives in the y and z directions and three unknown 
G (or E) components to determine at each grid point. As mentioned in 
Sec. |7.2.2[ it turns out that accurate center-difference approximations for 
the V X V X operator in three dimensions are better suited to a "staggered" 
grid, called a Yee grid |2,49 , in which different field components are dis- 
cretized at points slightly offset in space: e.g., Ex{[nx + h]Ax^nyAy^nzAz)^ 
Ey{nxAx^ [ny^^]Ay^ UzAz),^ and Ez{nxAx^ ^yAy^ [nz^\]Az) for the E field 
components. Note that any arbitrary frequency dependence of e is trivial to 
include, because in frequency domain one is solving each ^ separately, and a 
perfect electric conductor is simply the e{iC) -^ ^^ limit. 

One must, of course, somehow truncate the computational domain to a 
finite region of space in order to obtain a finite number N of degrees of free- 
dom. There are many reasonable ways to do this because Casimir interactions 
are rapidly decaying in space (force ~ l/a^+^ or faster with distance a in d 
dimensions, at least for zero temperature). One could simply terminate the 
domain with Dirichlet or periodic boundary conditions, for example, and if 
the cell boundaries are far enough away from the objects of interest then the 
boundary effects will be negligible (quite different from classical EM problems 
at real uj\) |12 . In classical EM, one commonly uses the more sophisticated 
approach of a perfectly matched layer (PML), an artificial refiectionless 
absorbing material placed adjacent to the boundaries of the computational 
domain to eliminate outgoing waves 1 2 . Mathematically, a PML in a direction 
X is equivalent to a complex "coordinate stretchin g" -^ -^ (1 -\-ia/uj)~ ^ 
for an artificial PML "conductivity" a{x) > (2l [9^[lQQ] , where the l/cu 



factor is introduced to give an equal attenuation rate at all frequencies. How- 
ever, at imaginary frequencies cj = i^, a PML therefore results simply in a 
real coordinate stretching ^ -^ (l + cr/^)~ ^: in a PDE with decaying, 
non-oscillatory solutions (such as the imaginary-cj Maxwell equations), it is 
well known that a reasonable approach to truncating infinite domains is to 
perform a (real) coordinate transformation that compresses space far away 
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where the solution is smah [53]. A convenience of Maxweh's equations is that 
any coordinate transformation (real or complex) can be converted merely into 
a change of e and ji flOl], so any PML can be expressed simply as a change 



of materials while keeping the same PDE and discretization |99, 100 . 

Such a center-difference scheme is nominally second-order accurate, with 
discretization errors that vanish as 0{Ax^) [2J|50]. One can also construct 
higher-order difference approximations (based on more grid points per differ- 
ence). As a practical matter, however, the accuracy is limited instead by the 
treatment of material interfaces where e changes discontinuously. If no spe- 
cial allowance is made for these interfaces, the method still converges, but its 



convergence rate is reduced by the discontinuity to 0{Ax) [54,55,102 103 



(unless one has E polarization completely parallel to all interfaces so that 
there is no field discontinuity). There are various schemes to restore second- 
order (or higher) accuracy by employing specialized FD equations at the 



interfaces [54} 142|, but an especially simple scheme involves unmodified FD 



equations with modified materials: it turns out that, if the discontinuous £ is 
smoothed in a particular way (to avoid introducing first-order errors by the 
smoothing itself), then second-order accuracy can be restored I55iil02i|103j P] 
Given the FD equations, one must must then choose a solution technique 
to solve the resulting linear equations Ax = b, where x is the Green's func- 
tion (or field), b is the delta- function (or current) right-hand side, and A is 
the discretized V x /i~^V x +^^£ operator. Note that A is a real-symmetric 



positive-definite matrix at imaginary frequencies, as discussed in Sec. |7.5.L2 
Because A is sparse [only 0{N) nonzero entries], one can utilize a sparse- 
direct Cholesky factorization A = R^R {R is upper-triangular) |66 (for 
which many software packages are available |68 ). Given this factorization, 
any right-hand side can be solved quickly by backsubstitution, so one can 
quickly sum the energy density over all grid points (essentially computing 
the trace of A~^) to find the Casimir energy, or alternatively sum the stress 
tensor over a stress-integration surface to find the force. Precisely such a 
sparse-direct FD method for the Casimir energy was suggested by Pasquali 
and Maggs |23 , albeit derived by a path-integral logdet expression that is 



mathematically equivalent to summing the energy density (7.8) |87 . The al- 
ternative is an iterative technique, and in this case A's Hermitian definiteness 
means that an ideal Krylov method, the conjugate-gradient method |65,67 
can be employed |12 . The conjugate-gradient method requires 0{N) stor- 
age and time per iteration, and in the absence of preconditioning requires a 
number of iterations in d dimensions proportional to the diameter 0{N^^^) 
of the grid for each right-hand side [104j . The stress-tensor approach reduces 



the number of right-hand sides to be solved compared to energy-density in- 



^ Even if the £ discontinuities are dealt with in this way, however, one may still fail 
to obtain second-order accuracy if the geometry contains sharp corners, which limit the 
accuracy to O(AxP) for some 1 < p < 2 102 . This is an instance of Darboux's principle: 
the convergence rate of a numerical method is generally limited by the strongest singularity 
in the solution that has not been explicitly compensated for [53]. 
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tegration: one only needs to evaluate the Green's function for sources on a 
stress-integration surface, which has 0{N^~) points in d dimensions. This 
gives a total time complexity of 0{N) • ©(A/"^/^) • 0(A/''^) = 0{N'^) for an 
unpreconditioned iterative method; an ideal multigrid preconditioner can in 
principle reduce the number of iterations to 0(1) |4 |lQ5| (when N is increased 
by improving spatial resolution), yielding an 0{N'^~^) time complexity. Sub- 
stantial further improvements are obtained by realizing that one does not, in 
fact, need to sum over every point on the stress-integration surface, instead 
switching to a different spatial integration scheme described in Sec. |7.5.6[ 



7.5.5 Boundary- element methods (BEMs) 

In some sense, a volume discretization such as an FD method is too general: 
in most physical situations, the medium is piecewise-constant, and one might 
want to take advantage of this fact. In particular, for the basic problem of 
finding the field in response to a current source at a given frequency, one can 
instead use a surface integral-equation approach: the unknowns are surface 
currents on the interfaces between homogeneous materials, and one solves for 
the surface currents so that the total field (source + surface currents) satisfies 
the appropriate boundary conditions at the interfaces [l]|3][7]. For example, in 
the case of a perfect electric conductor, the surface-current unknowns can be 
the physical electric currents J at the interface, and the boundary condition 
is that of vanishing tangential E field {^ In the case of permeable media e 
and /i, the physical (bound) currents are volumetric within the medium [e.g., 
the electric bound current is J = — icj(£ — £o)E], not surface currents [4l] . 
However, it turns out that one can introduce fictitious surface electric and 
magnetic currents at all interfaces to provide enough degrees of freedom to 
satisfy the boundary condition of continuous tangential E and H, and thus 
to fully solve Maxwell's equations. The application of this equivalence princi- 
pl^^to obtain surface integral equations for BEM is known as the PMCHW 



approach (Poggio, Miller, Chang, Harrington, and Wu) |11Q112| . In either 



case, one has surface (electric and/or magnetic) currents Jg, plus an external 



current source J [e.g., the right-hand side of (7.4)], so one can express the E or 
H field at any point x as a convolution of J + Jg with the analytically known 
Green's function Go(x — x') of the corresponding homogeneous medium at 
X. In BEM, one expresses J^, in turn, as a sum of localized basis functions h^ 



^^ This is known as an electric-field integral equation (EFIE); one can also express the 
equations for perfect conductors in terms of boundary conditions enforced on magnetic 
fields (MFIE) or some linear combination of the two (CFIE), and the most effective 
formulation is still a matter of debate |90|. 

^^ The idea of solving scattering problems by introducing fictitous boundary currents had 
its origins |106Hl09] many years before its application to BEM by Harrington [110] and 
subsequent refinements. 
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Fig. 7.3 Example triangular mesh of the surfaces of two objects for a BEM solver 
Associated with each edge k is an "RWG" basis function b^ |63 , schematically represented 
in the inset, which vanishes outside the adjacent two triangles. 



associated with some discrete mesh approximation of the surface. For exam- 
ple, Fig. |7.3| depicts a standard triangular-type mesh of two objects, where 
there is a localized basis function h^ (inset) associated with each edge of this 
mesh such that b/^ is nonzero only on the adjacent two triangles |63 ; this is 
the RWG basis mentioned in Sec. |7.2.2[ Abstractly, the resulting equations 
for the fields could then be written in the following form: 

N 

field(x) = Go * (J + J.) = Go * J + ^ Go * b^c^, (7.16) 

k=l 

where Go* denotes convolution with the (dyadic) analytical Green's function 
of the homogeneous medium at x, and c^ are the unknown coefficients of each 
basis function}^ (More generally. Go * J could be replaced by any arbitrary 
incident field, regardless of how it is created.) In a Galerkin method (see 



Sec. 7.2.2), one obtains N equations for the TV unknowns Ck by taking the 



inner product of both sides of this equation (substituted into the appropriate 
boundary condition) with the same basis functions hj (since they work just as 
well as a basis for the tangential field as for the tangential surface currents). 
This ultimately results in a set of linear equations Ac = d, where the matrix 
A multiplying the unknown coefficients Ck is given by 



Technically, only currents from surfaces bordering the medium of x contribute to this 
sum. 
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^jk= ffb,(x)-Go(x-xO-bfc(xO^'xdV. (7.17) 

[For the case of a perfect conductor with a vanishing tangential E, the right- 
hand-side d is given by dj = — (bj,Go * J) = — ^bj(x) • Go(x — x') • 
3(yi') dpyidpyi' .] One then solves the linear system for the unknown coeffi- 
cients c/e, and hence for the unknown surface currents Jg. Implementing this 
technique is nontrivial because the Ajj. integrands (7.17) are singular for 



j = k OT for j adjacent to /c, necessitating specialized quadrature techniques 
for a given form of Go |6Q|[61] , but substantial guidance from the past decades 
of literature on the subject is available. 

Given these currents, one can then evaluate the electric or magnetic field at 



any point x, not just on the surface, by evaluating eq. (7.16) at that point. In 



particular, one can evaluate the field correlation functions via the fluctuation- 



dissipation theorem (7.3): {Ej{-K)Ek{'x.))^ is given in terms of the electric field 
in the j direction at x from a delta-function current J in the k direction at 
X. Of course, as noted previously, this is infinite because the Go * J term (the 
field from the delta function) blows up at x, but in the Casimir case one is only 
interested in the change of the correlation functions due to the geometry — so, 
one can use the standard trick [81 1 of subtracting the vacuum contribution 
Go * J and only computing the surface-current contribution Go * Js to the 
field at x. In this way, one can compute the stress tensor, the energy density, 
and so on, as desired. 

As explained in Sec. |7.4[ the integral of contributions over all frequencies is 
best performed at imaginary frequencies, so all of the above must use uj = i^. 
This only has the effect of Wick-rotating the homogeneous-medium dyadic 
Green's function Go to the ^ e~^^^~^ l/|x — x'| imaginary-frequency Green's 
function. This makes the problem easier, in principle. First, the exponen- 
tial decay cuts off long-range interactions, making fast-solver techniques (see 



Sec. 7.2.3) potentially even more effective. Second, the matrix A is now real- 
symmetric and positive-definite, which allows the use of more efficient linear 
solvers as noted previously. Fortunately, the l/|x — x'| singularity of Go is 
the same at real and imaginary frequencies, allowing existing techniques for 



the integration of (7.17) to be leveraged. 

At first glance, this approach seems most straightforwardly applicable to 
the stress- tensor technique, as suggested in Ref. |12 : one uses the BEM solu- 
tion to evaluate the mean stress tensor (T) on any integration surface around 
a body, integrating via some quadrature technique to obtain the force. If one 
uses a dense-direct solver (when N is not too big), the Cholesky factoriza- 
tion of A can be computed once for a given ^ and then many right-hand 
sides can be solved quickly via backsubstitution |65 in order to integrate 
(T) over the stress-integration surface. Precisely such a dense-direct BEM 
stress-tensor method was recently demonstrated to compute Casimir forces 
in two and three dimensions y^|20|. As described in Sec. |7.2.3[ fast-solver 
techniques can be applied to multiply A by a vector in O(A^logA^) time 
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with 0{N) storage; given a good preconditioner, this imphes that an itera- 
tive method such as conjugate-gradient (apphcable since A is real-symmetric 
positive-definite) could find (T) at a single x and ^ in 0{N log N) time. The 
remaining question is the number of points required for the surface integral 
of (T), which depends on why one is increasing TV: either to increase accuracy 
for a fixed geometry or to increase the complexity of the geometry for a fixed 
accuracy. In the former case, the smoothness of (T) in x means that exponen- 
tially convergent quadrature techniques are applicable, which converge much 
faster than the (polynomial) BEM basis for the surface currents, so that ul- 
timately the number of stress-quadrature pointq^ should be independent of 
N and the overall complexity becomes O(A^logA^). In the latter case, for a 
fixed accuracy and increasingly complex geometry (or smaller feature sizes), 
it appears likely that the number of stress-quadrature points will increase 
with A/", but detailed studies of this scaling are not yet available. 

It turns out that this BEM approach is closely related to the BEM path- 
integral approach described in Sec. |7.6.3| Both approaches end up solving 



linear equations with exactly the same matrix A of eq. (7.17), with the same 
degrees of freedom. The path-integral approach shows, however, that this 
same matrix can be applied to compute the Casimir interaction energy as well 
as the force, with comparable computational cost for dense solvers. Moreover, 
as explained below, expressing the force in terms of the derivative of the path- 
integral energy results in a trace expression that is conceptually equivalent 
to integrating a stress tensor over the surface of an object, where the number 
of "quadrature points" is now exactly equal to A^. An unanswered question, 
at this point, is whether a fast solver can be more efficiently (or more easily) 
exploited in the stress-tensor approach or in the path-integral approach. 



7.5.4 Other possibilities: FEM and spectral methods 

There are of course, many other frequency- domain techniques from classical 
EM that could potentially be used to solve for the Green's function and 
hence the energy /force density. For example, one could use spectral integral- 
equation methods, such as multipole expansions for spheres and cylinders [5], 
to compute responses to currents, although the advantages of this approach 
compared to the spectral path-integral approach in Sec. |7.6.2| are unclear. 



^^ Numeric integration {quadrature) approximates an integral f f{x)dx by a sum 
y^ ^- f(xi)wi over quadrature points Xi with weights Wi. There are many techniques for the 
selection of these points and weights, and in general one can obtain an error that decreases 
exponentially fast with the number of points for analytic integrands |53|83[84[|113] . Multi- 
dimensional quadrature, sometimes called cubature, should be used to integrate the stress 
tensor over a 2d surface, and numerous schemes have been developed for low-dimensional 
cubature ^114|115] (including methods that adaptively place more quadrature points where 
they are most needed |116) ). For spherical integration surfaces (or surfaces that can be 
smoothly mapped to spheres), specialized methods are available (11 711 1181. 
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One can also solve the PDE formulation of the Green's function (7.10) using a 
finite-element (FEM) approach with some general mesh; in principle, existing 
FEM techniques from classical EM flJISpllT] are straightforwardly applicable. 
One subtlety that arises in FEM methods with a nonuniform resolution is the 
regular izat ion, however |12 . In principle, as mentioned above, one needs to 
subtract the vacuum Green's function contribution from the field correlation 
functions in order to get a physical result [since the vacuum Green's function 
G(x, x') diverges as x' -^ x, although the divergence is cut off by the the finite 
mesh resolution]. With a uniform mesh, this vacuum contribution is the same 
everywhere in the mesh and hence automatically integrates to zero in the force 
(when the stress tensor is integrated over a closed surface or the energy is 
differentiated). For a nonuniform mesh, however, the vacuum contribution 
varies at different points in space with different resolution, so some "manual" 
regularization seems to be required (e.g., subtracting a calculation with the 
same mesh but removing the objects). These possibilities currently remain 
to be explored for Casimir physics. 



7.5.5 Finite- difference time-domain (FDTD) methods 

Casimir effects are fundamentally broad-bandwidth, integrating contribu- 
tions of fluctuations at all frequencies (real or imaginary), although the 
imaginary-frequency response is dominated by a limited range of imaginary 
frequencies. In classical EM, when a broad-bandwidth response is desired, 
such as a transmission or reflection spectrum from some structure, there is 
a well-known alternative to computing the contributions at each frequency 
separately — instead, one can simulate the same problem in time^ Fourier- 
transforming the response to a short pulse excitation in order to obtain the 
broad-bandwidth response in a single time- domain simulation [6, 119|. The 



same ideas are applicable to the Casimir problem, with a few twists, yield- 
ing a practical method [l3J[27] that allows Casimir calculations to exploit 
off-the-shelf time-domain solvers implementing the standard finite- difference 
time-domain (FDTD) method [2]. There are two key components of this ap- 
proach |13|: first, converting the frequency integral to a time integral and, 
second, finding a time-domain equivalent of the complex- frequency idea from 

Sec. [721 

As reviewed above, the mean fluctuations in the fields, such as (£^^(x))c^, 
can be expressed in terms of the fields at x from a frequency-cj current at 
X. If, instead of a frequency-o; current, one uses a current with 5{t) time 



dependence, it follows by linearity of (7.4) that the Fourier transform of the 
resulting fields must yield exactly the same {E'^{^))^. Roughly, the procedure 
could be expressed as follows: First, we compute some function r{t) of the 
time-domain fields from a sequence of simulations with 5(t) sources, e.g. where 
r{t) is the result of spatially integrating the fields making up the mean 
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stress tensor (T(x)) [noting that each point x involves several separate S{t)- 
response simulations]. Second, we Fourier transform r{t) to obtain r{uj). 
Third, we obtain the force (or energy, etcetera) by integrating f r{uj)g{uj)duj 
with appropriate frequency- weighting factor g{uj) (which may come from the 
frequency dependence of e in (T), a Jacobian factor from below, etcetera). 
At this point, however, it is clear that the Fourier transform of F was entirely 
unnecessary: because of the unitarity of the Fourier transform (the Plancherel 
theorem), J r{uj)g{uj)duj = J r{t)g{—t)dt. That is, we can compute the force 
(or energy, etcetera) by starting with S{t) sources and simply integrating the 
response r{t) in time (accumulated as the simulation progresses) multiplied 
by some (precomputed, geometry-independent) kernel g{t) (which depends 
on temperature if the coth factor is included for T > 0) . The details of this 
process, for the case of the stress tensor, are described in Refs. [l3J[27l. 

Although it turns out to be possible to carry out this time-integration 
process as-is, we again find that a transformation into the complex- frequency 
plane is desirable for practical computation (here, to reduce the required 
simulation time) [iS]. Transforming the frequency in a ^ime-domain method, 
however, requires an indirect approach. The central observation is that, in the 



equation (7.4) for the electric-field Green's function G^ , the frequency only 
appears explicitly in the uj'^e term, together with e. So, any transformation of 
uj can equivalently be viewed as a transformation of e. In particular suppose 
that we wish to make some transformation cu -^ co'(^) to obtain an u in 
the upper-half complex plane, where ^ is a real parameter (e.g. u = i(_ for 
a Wick rotation). Equivalently, we can view this as a calculation at a real 
frequency ^ for a transformed complex material: uj^£{uj, x) -^ £,'^£c{£,^ x) where 



the transformed material is 13, 120 






e,(^,x) = ^eMC),x). (7.18) 



For example, a Wick rotation u ^ i(_ is equivalent to operating at a real 
frequency ^ with a material e{uj) -^ —e{i(,). However, at this point we run 
into a problem: multiplying s hy —1 yields exponentially growing solutions 
at negative frequencies [l3J[l20|, and this will inevitably lead to exponen- 
tial blowup in a time-domain simulation (which cannot avoid exciting neg- 
ative frequencies, if only from roundoff noise). In order to obtain a useful 
time-domain simulation, we must choose a contour a;(^) that yields a causal, 
dissipative material Sc^ and one such choice is uj{^) = ^y^T^Ho/^ for any 
constant a > 13, 12Q|. This yields £c = {I -\- icr/uj)£, where the ia/uj term 



behaves exactly like an artificial conductivity added everywhere in space. In 



the frequency-domain picture, we would say from Sec. 7.4 that this uj{^) 
contour will improve the computation by moving away from the real-cj axis, 
transforming the frequency integrand into something exponentially decaying 
and less oscillatory. In the time-domain picture, the a term adds a dissipation 
everywhere in space that causes r{t) to decay exponentially in time, allow- 
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ing us to truncate the simulation after a short time. As long as we include 
the appropriate Jacobian factor ^ in our frequency integral, absorbing it 
into g{t), we will obtain the same result in a much shorter time. The com- 
putational details of this transformation are described in Refs. |13,27 . More 
generally, this equivalence between the Casimir force and a relatively narrow- 
bandwidth real-frequency response of a dissipative system potentially opens 
other avenues for the understanding of Casimir physics |120 . 

The end result is a computational method for the Casimir force in which 
one takes an off-the-shelf time-domain solver (real time/frequency) , adds an 
artificial conductivity a everywhere, and then accumulates the response r{t) 
to short pulses multiplied by a precomputed (geometry independent) kernel 
g{t). The most common time-domain simulation technique in classical EM is 
the FDTD method |2 . Essentially, FDTD works by taking the same spatial 
Yee discretization as in the FDFD method above, and then also discretizing 
time with some time step At. The fields are then marched through time in 
steps of At^ where each time step requires 0{N) work for N spatial grid 
points. Because the complex-cj contour is implemented entirely as a choice 
of materials Sc, existing FDTD software can be used without modification 
to compute Casimir forces, and one can exploit powerful existing software 
implementing parallel calculations, various dimensionalities and symmetries, 
general dispersive and anisotropic materials, PML absorbing boundaries, and 
techniques for accurate handling of discontinuous materials. One such FDTD 
package is available as free/open-source software from our group |119| , and 
we have included built-in facilities to compute Casimir forces fl211. 



7.5.6 Accelerating FD convergence 

Finally, we should mention a few techniques that accelerate the convergence 
and reduce the computational cost of the finite-difference approaches. These 
techniques are not necessary for convergence, but they are simple to imple- 
ment and provide significant efficiency benefits. 

The simplest technique is extrapolation in Ax: since the convergence rate 
of the error with the spatial resolution Ax is generally known a priori, one 
can fit the results computed at two or more resolutions in order to extrapolate 
to Ax -^ 0. The generalization of this approach is known as Richardson ex- 
trapolation [91], and it can essential increases the convergence order cheaply, 
e.g., improving 0{Ax) to 0{Ax^) |123 . 

Second, suppose one is computing the force between two objects A and 
B surrounded by a homogeneous medium. If one of the objects, say B, is 
removed, then (in principle) there should be no net remaining force on A. 
However, because of discretization asymmetry, a computation with A alone 
will sometimes still give a small net force, which converges to zero as Ax -^ 0. 
If this "error" force is subtracted from the A-B force calculation, it turns out 
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that the net error is reduced. More generahy, the error is greatly reduced if 
one computes the A-B force and then subtracts the "error" forces for A alone 
and for B alone, tripling the number of computations but greatly reducing 
the resolution that is required for an accurate result |12 . 

Third, when integrating the stress tensor (T(x))^^ over x to obtain the 



net force (7.12), the most straightforward technique in FD is to simply sum 
over all the grid points on the integration surface — recall that each point x 
requires a linear solve (a different right-hand side) in frequency domain, or 
alternatively a separate time-domain simulation (a separate current pulse). 
This is wasteful, however, because (T(x))^^ is conceptually smoothly varying 
in space — if one could evaluate it at arbitrary points x (as is possible in 
the BEM approach), an exponentially convergent quadrature scheme could 
be exploited to obtain an accurate integral with just a few x's. This is not 
directly possible in an FD method, but one can employ a related approach. 
If the integration surface is a box aligned with the grid, one can expand the 
fields on each side of the box in a cosine series (a discrete cosine transform, or 
DCT, since space is discrete) — this generally converges rapidly, so only a small 
number terms from each side are required for an accurate integration. But 
instead of putting in point sources, obtaining the responses, and expanding 
the response in a cosine series, it is equivalent (by linearity) to put in cosine 
sources directly instead of point sources. [Mathematically, we are exploiting 
the fact that a delta function can be expanded in any orthonormal basis 6n(x) 
over the surface, such as a cosine series, via: (5(x — x') = Xln ^^(^O^n(x). 



Substituting this into the right-hand side of (7.10), each ^^(x) acts like a 
current source and bn{'x.') scales the result, which is eventually integrated 
over x^] The details of this process and its convergence rate are described 
in Ref. [27], but the consequence is that many fewer linear systems (fewer 
right-hand sides) need be solved (either in frequency or time domain) than 
if one solved for the stress tensor at each point individually. 



7.6 Path integrals and scattering matrices 

Another formulation of Casimir interactions is to use a derivation based on 
path integrals. Although the path-integral derivation itself is a bit unusual 
from the perspective of classical EM, and there are several slightly different 
variations on this idea in the literature, the end result is straightforward: 
Casimir energies and forces are expressed in terms of log determinants and 
traces of classical scattering matrices [10,14-18,21,22,24-26], or similarly 



the interaction matrices (7.17) that arise in BEM formulations [9]. Here, we 



omit the details of the derivations and focus mainly on the common case of 
piecewise- homogeneous materials, emphasizing the relationship of the result- 
ing method to surface-integral equations from classical EM via the approach 
in Ref. 19. 
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Path integrals relate the Casimir interaction energy /7 of a given config- 
uration to a functional integral over all possible vector-potential fields A. 
Assuming piecewise- homogeneous materials, the constraint that the fields 
in this path integral must satisfy the appropriate boundary conditions can 
be expressed in terms of auxiliary fields J at the interfaces (a sort of La- 
grange multiplier) |124 F^At this point, the original functional integral over 
A can be performed analytically, resulting in an energy expression involving 
a functional integral Z{^) over only the auxiliary fields J at each imaginary 
frequency ^ , of the form (at zero temperature) : 

^.-|fl„g*.£|^. (7.19, 

^r^\ ^ / pjg-^//(i^x//d2x'j(x)-G5(x-x')-J(x')^ (7.20) 

Here, Z^o denotes Z when the objects are at infinite separation (non- 
interacting), regularizing U to just the (finite) interaction energy (See also 
the chapter by S. J. Rahi et al. in this volume for additional discussion of 
path integrals and Casimir interactions.) In the case of perfect electric con- 
ductors in vacuum, J can be interpreted as a surface current on each con- 
ductor (enforcing the vanishing tangential E field), and G^ is the vacuum 
Green's function in the medium outside the conductors |9 . For permeable 
media (finite e and /i), it turns out that a formulation closely related to the 



standard PMCHW integral-equation model (see Sec.|7.5.3) can be obtained: 



J represents fictitious surface electric and magnetic currents on each interface 
(derived from the continuity of the tangential E and H fields), with G^ again 
being a homogeneous Green's function (with one Z factor for each contiguous 
homogeneous region) [125J . Alternatively, because there is a direct correspon- 
dence between surface currents and the outgoing/scattered fields from a given 
interface, "currents" J can be replaced by scattered fields, again related at dif- 
ferent points X and x' by the Green's function of the homogeneous medium; 
this is typically derived directly from a T-matrix formalism p!4, 22, 25, 26], 



Here, we will focus on the surface-current viewpoint, which is more common 
in the classical-EM integral-equation community. 



The path integral (7.20) is somewhat exotic in classical EM, but it quickly 
reduces to a manageable expression once an approximate (finite) basis h^ is 
chosen for the currents J. Expanding in this basis, J ~ ^c/eb/e(x) and the 
functional integral VJ is replaced by an ordinary integral over the basis co- 



efficients dci • • • dcN. Equation (7.20) is then a Gaussian integral that can be 



performed analytically to obtain Z{^) = #/^/detA{£^) for a proportionality 
constant # |9 , where Ajk = fhj - G^ - hk is essentially the same as the 



^^ Alternatively, the path integral can be performed directly in A, resulting in an expression 
equivalent to the sum over energy density in Sec. |7.5| |87 and which in an FD discretization 
reduces in the same way to repeated solution of the Green's-function diagonal at every point 
in space |23|. 
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BEM matrix (7.17), albeit here in an arbitrary basis. In the logdet of (7.19) 



proportionahty constants and exponents cancel, leaving: 

he 



^ '2. 



/ logdet [Aoo(0"'^(0]^e (7.21) 



Just as in Sec. 7.17[ the use of a real-symmetric positive-definite homogeneous 



Green's function G^ at imaginary frequencies means that A(^) is also real- 
symmetric and positive-definite, ensuring positive real eigenvalues and hence 
a real log det. Several further simplifications are possible, even before choosing 
a particular basis. For example, let p be the position of some object for which 
the force F is desired. The components Fi of the force (in direction pi) can 
then be expressed directly as a trace [9 , 126 : 



/ tr(A-^^)de (7.22) 



dpi 27r io V ^P 



Equivalently, this trace is the sum of eigenvalues A of the generalized eigen- 
problem ^v = AAv; again, these A are real because A is real-symmetric 
positive-definite and dA/dpi is real-symmetric. (If dense-direct solvers are 
used, computing A~^^ via Cholesky factorization is much more efficient 
than computing eigenvalues, however |65 .) The matrix A can be further 
block-decomposed in the usual case where one is computing the interactions 
among two or more disjoint objects (with disjoint surface currents J). For 
example, suppose that one has two objects 1 and 2, in which case one can 
write 

where An and A22 couple currents on each object to other currents on the 
same object, and ^412 and AJ2 = ^21 couple currents on object 1 to object 2 
and vice versa. In the limit of infinite separation for Aoo, one obtains ^412 -^ 



while All and ^22 are unchanged, and one can simplify the logdet in (7.21) 
to 

logdet [A^(0-'A(0] = logdet [/ - ^"2^^^ ^-^^21] . (7.24) 

Computationally, only A12 depends on the relative positions of the objects, 
and this simplification immediately allows several computations to be re-used 
if the energy or force is computed for multiple relative positions. 



7.6.1 Monte- Carlo path integration 

Before we continue, it should be noted that there also exists a fundamentally 
different approach for evaluating a path-integral Casimir formulation. Instead 
of reducing the problem to surface/scattering unknowns and analytically in- 
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tegrating Z to obtain a matrix log det expression, it is possible to retain the 
original path-integral expression, in terms of a functional integral over vector 
potentials A in the volume, and perform this functional integral numerically 
via Monte-Carlo methods [127, 1281. This reduces to a Monte-Carlo integra- 
tion of an action over all possible closed-loop paths ("worldlines") , discretized 
into some number of points per path. Because this technique is so different 
from typical classical EM computations, it is difficult to directly compare 
with the other approaches in this review. Evaluating its computational re- 
quirements involves a statistical analysis of the scaling of the necessary num- 
ber of paths and number of points per path with the desired accuracy and the 
complexity of the geometry |12 , which is not currently available. A difficulty 
with this technique is that it has currently only been formulated for scalar 
fields with Dirichlet boundary conditions, not for the true Casimir force of 
vector electromagnetism. 



1.6.2 Spectral methods 

One choice of basis functions hk for the path-integral expressions above is 
a spectral basis, and (mirroring the history of integral equations in classi- 
cal EM) this was the first approach applied in the Casimir problem. With 
cylindrical objects, for example, the natural spectral basis is a Fourier-series 
g*m0 ^^ ^Yie angular direction ^. For planar surfaces the natural choice is a 
Fourier transform, for spheres it is spherical harmonics Yim (or their vec- 
tor equivalents |41 ), and for spheroids there are spheroidal harmonics [129|. 



Equivalently, instead of thinking of surface currents expanded in a Fourier-like 
basis, one can think of the scattered fields from each object expanded in the 
corresponding Fourier-like basis (e.g. plane, cylindrical, or spherical waves), 
in which case A relates the incoming to outgoing/scattered waves for each 
object; this has been called a "scattering-matrix" or "T- matrix" method and 
is the source of many pioneering results for Casimir interactions of non-planar 
geometries [14,22,25,26 . Even for nonspherical/spheroidal objects, one can 



expand the scattered waves in vector spherical harmonics |22| , and a variety 
of numerical techniques have been developed to relate a spherical-harmonic 
basis to the boundary conditions on nonspherical surfaces |58 . These spec- 
tral scattering methods have their roots in many classical techniques for EM 
scattering problems [56 , 122| (See also the chapters by S. J. Rahi et al. and 



A. Lambrecht et al. in this volume for additional discussions of scattering 
techniques and Casimir interactions.) Here, we will use the surface-current 
viewpoint rather than the equivalent scattered- wave viewpoint. 



Many simplifications occur in the interaction matrix A of (7.23) for geome- 
tries with highly symmetrical objects and a corresponding spectral basis [22] . 
Consider, for example, the case of spherical objects, with surface currents ex- 
pressed in a vector spherical-harmonic basis (spherical harmonics for two po- 
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larizations i22i). In the interaction matrix Ajk = ff bj(x) • G(x — x') • b/i;(x'), 
the convolution f G(x — x') •b/e(x') of a Green's function G with h^ is known 
analyticahy: it is just the outgoing spherical wave produced by a spherical- 
harmonic current. If bj is another spherical-harmonic current on the same 
sphere, then the orthogonality of the spherical harmonics means that the x 
integral of hj (x) against the spherical wave is zero unless j = k. Thus, the self- 



interaction blocks All and A22 of (7.23), with an appropriate normalization 



are simply identity matrices. The ^412 entries are given by the coupling of a 
spherical wave from b^ on sphere 2 with a spherical-harmonic basis function 
hj on sphere 1, but again this integral can be expressed analytically, albeit 
as an infinite series: the spherical wave from sphere 2 can be re-expressed in 
the basis of spherical waves centered on sphere 1 via known translation iden- 
tities of spherical waves, and as a result ^412 takes the form of a "translation 
matrix" |22|. Furthermore, if there are only two spheres in the problem, then 
their spherical harmonics can be expressed with respect to a common z axis 
passing through the centers of the spheres, and a Y^m on sphere 1 will only 
couple with a Y^/m' on sphere 2 if m = m', greatly reducing the number of 
nonzero matrix elements. Related identities are available for coupling cylin- 
drical waves around different origins, expanding spherical/cylindrical waves 
in terms of planewaves for coupling to planar surfaces, and so on |22|. 

As was noted in Sec. |7.2.2[ such a spectral basis can converge exponen- 
tially fast if there are no singularities (e.g. corners) that were not accounted 
for analytically, and the method can even lend itself to analytical study. Espe- 
cially for cylinders and spheres, the method is simple to implement and allows 
rapid exploration of many configurations; the corresponding classical "multi- 
pole methods" are common in classical EM for cases where such shapes are of 
particular interest ^. On the other hand, as the objects become less and less 
similar to the "natural" shape for a given basis (e.g. less spherical for spher- 
ical harmonics), especially objects with corners or cusps, the spectral basis 
converges more slowly f58^. Even for the interaction between two spheres or 
a sphere and a plate, as the two surfaces approach one another the multipole 
expansion will converge more slowly [25[|13Q[p^ — conceptually, a spherical- 



harmonic basis has uniform angular resolution all over the sphere, whereas 
for two near-touching surfaces one would rather have more resolution in the 
regions where the surfaces are close (e.g. by using a nonuniform BEM mesh). 
This exponential convergence of a spectral (spherical harmonic |22 ) Casimir 
calculation is depicted in Fig. |7.4| for the case of the Casimir interaction 
energy U between two gold spheres of radius R = 1 /im, for various surface- 
to-surface separations a. The error AU/U decreases exponentially with the 
maximum spherical-harmonic order i [corresponding to N = 4:i{i-\- 2) degrees 
of freedom for two spheres], but the exponential rate slows as a/R decreases. 
(On the other hand, for small a/R a perturbative expansion or extrapolation 
may become applicable |13Q|.) 
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Fig. 7.4 Fractional error AU/U in the Casimir interaction energy U between two gold 
spheres of radius R = 1 /im, for various surface-surface separations a, as a function 
of the maximum spherical-harmonic order i of the spectral path-integral (scattering- 
matrix/multipole) method. The error converges exponentially with i, but the exponential 
rate slows as a/R shrinks. (Calculations thanks to A. Rodriguez.) 



7.6.3 Boundary- element methods (BEMs) 



In a BEM, one meshes the interfaces, say into triangles, and uses a set of 
localized basis functions b/^ as discussed in Sec. |7.5.3[ In this case, the in- 
teraction matrix A that arises in the path-integral formulation is exactly the 
same as the interaction matrix that arises in classical BEM methods (albeit 
at an imaginary frequency), and is the same as the matrix A that arises in a 



BEM stress-tensor approach as described in Sec. |7.5.3| The main difference, 
compared to the stress-tensor approach, lies in how one uses the matrix A\ 
instead of solving a sequence of linear equations to find the mean stress tensor 
(T) at various points on a surface around an object, one computes logdet A 

or tr A~^ ^ to obtain the energy (7.21) or force (7.22). We have demon- 
strated this approach for several three-dimensional geometries, such as the 
crossed capsules of Fig. 7.3 |9 . 



If one is using dense-matrix techniques, the advantage of this approach 
over the stress-tensor technique seems clear |9 : it avoids the complication 
of picking a stress-integration surface and an appropriate surface-integration 
technique, and allows the size of the linear system to be easily reduced via 



blocking as in eq. (7.24). The situation is less clear as one moves to larger and 



larger problems, in which dense-matrix solvers become impractical and one 
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requires an iterative method. In that case, computing tr A~^^ straight- 
forwardly requires N hnear systems to be solved; if each linear system can 



be solved in 0(A^log A^) time with a fast solver (as discussed in Sees. 7.2.3 
and 7.5.3), then the overall complexity is 0{N'^ log N) [with 0{N) storage]. 



whereas it is possible that the stress-tensor surface integral may require fewer 
than N solves. On the other hand, there may be more efficient ways to com- 
pute the trace (or log det) via low-rank approximations: for example, if the 
trace (or log det) is dominated by a small number of extremal eigenvalues, 
then these eigenvalues can be computed by an iterative method |68 with the 
equivalent of <C A/" linear solves. The real-symmetric property of A, as usual, 
means that the most favorable iterative methods can be employed, such as 
a Lanczos or Rayleigh-quotient method |68 . Another possibility might be 
sparse-direct solvers via a fast-multipole decomposition |73 . The most effi- 
cient use of a fast 0{N log N) BEM solver in Casimir problems, whether by 
stress-tensor or path-integral methods, remains an open question (and the 
answer may well be problem-dependent). 

In the BEM approach with localized basis functions, the tr A~^ ^ ex- 
pression for the force corresponds to a sum of a diagonal components for each 
surface element, and in the exact limit of infinite resolution (infinitesimal el- 
ements) this becomes an integral over the object surfaces. Expressing the 
force as a surface integral of a quantity related to Green's-function diagonals 
is obviously reminiscent of the stress-tensor integration from Sec. |7. 5.0} and 
it turns out that one can prove an exact equivalence using only vector cal- 
culus [125] . (At least one previous author has already shown the algebraic 
equivalence of the stress tensor and the derivative of the path-integral energy 
for forces between periodic plates [1311.) 



7.6.4 Hybrid BEM /spectral methods 

It is possible, and sometimes very useful, to employ a hybrid of the BEM and 
spectral techniques in the previous two sections. One can discretize a surface 
using boundary elements, and use this discretization to solve for the scatter- 
ing matrix A^^ of each object in a spectral basis such as spherical waves. That 
is, for any given incident spherical wave, the outgoing field can be computed 



with BEM via (7.16) and then decomposed into outgoing spherical waves to 
obtain one row/column of Akk at a time; alternatively, the multipole decom- 
position of the outgoing wave can be computed directly from the multipole 
moments of the excited surface currents J5 |41 . This approach appears to 
be especially attractive when one has complicated objects, for which a lo- 
calized BEM basis works well to express the boundary conditions, but the 
interactions are only to be computed at relatively large separations where 
the Casimir interaction is dominated by a few low-order multipole moments. 
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Fig. 7.5 Schematic problem for which eigenmode-expansion is well suited: the interaction 
between two corrugated surfaces, with period A. The Casimir problem reduces to comput- 
ing the reflection matrices Ri and R2 for each individual surface, in a planewave basis. 
Eigenmode expansion works by expanding the fleld in each cross-section (dashed lines) 
in the basis of eigenmodes of a z-invariant structure with that cross-section, and then 
matching boundary conditions whenever the cross-section changes. 



One can perform the BEM computation once per object and re-use the re- 
sulting scattering matrix many times via the analytical translation matrices, 
allowing one to efficiently compute interactions for many rearrangements of 
the same objects and/or for "dilute" media consisting of many copies of the 
same objects [132| . (Essentially, this could be viewed as a form of low-rank 
approximation of the BEM matrix, capturing the essential details relevant 
to moderate-range Casimir interactions in a much smaller matrix.) Such a 
hybrid approach is less attractive for closer separations, however, in which 
the increasing number of relevant multipole moments will eventually lead to 
an impractically large matrix to be computed. 



7.6.5 Eigenmode- expansion /RCW A methods 



Consider the case of the interaction between two corrugated surfaces depicted 



in Fig 7.5, separated in the z direction. From the scattering-matrix viewpoint, 
it is natural to consider scattering off of each object by plane waves. In this 
case, the self-interaction matrices A^^ and A22 can be re-expressed in terms 
of reflection matrices Ri and R2 for each surface, relating the amplitudes of 
incident waves at some plane (dashed line) above each surface to the reflected 
(specular and non-specular) planewave amplitudes. The matrices ^412 and 
A21 are replaced by a diagonal matrix D12 = ^21 ^^^^ relates the planewave 
amplitudes at the planes for objects 1 and 2, separated by a distance a — at 
real frequencies, this would be a phase factor, but at imaginary frequencies 
it is an exponential decay as discussed below. This results in the following 
expression for the Casimir interaction energy: 
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he r^^ 

U=— logdet [/ - R2Di2RiDi2] di. (7.25) 

Alternatively, instead of viewing it as a special case of the T-matrix/scattering- 
matrix idea [22], the same expression can be derived starting from an 
eigenmode-summation approach |18 . 

The problem then reduces to computing the scattering of an incident 
planewave off of a corrugated surface, with the scattered field decomposed 
into outgoing planewaves. For this problem, one could use any of the tools 
of computational EM (such as BEM, FD, and so on), but there is a no- 
table method that is often well-suited to the case of periodic surfaces, es- 
pecially periodic surfaces with piecewise-constant cross-sectionq^ (as in ob- 

or 



ject 2 of Fig. 7.5). This method is called eigenmode expansion [43f 45 
rigorous coupled-wave analysis (RCWA) |46,47|, or alternatively a cross- 
section method |48|. RCWA has a long history because it is closely tied 
to semi-analytical methods to study waveguides with slowly/weakly varying 
cross-sections [48" 133]. An analogous method was recently applied to Casimir 
problems |18 . In RCWA, one computes reflection and scattering matrices at 
a given frequency u along some direction z by expanding the fields at each z 
in the basis of the eigenmodes of the cross-section at that z (waves with z de- 
pendence e*^^ at a given cj, where f3 is called the propagation constant of the 
mode). Along regions of uniform cross-section, the z dependence e*^^ of each 
mode is known analytically and no computation is required (the mode am- 
plitudes are multiplied by a diagonal propagation matrix D). Regions of con- 
tinuously varying cross-section are approximated by breaking them up into 



a finite number of constant-cross-section layers (as in object 1 of Fig. 7.5). 
At any z where the cross-section changes, a change of basis is performed by 
matching boundary conditions (the xy components of the fields must be con- 
tinuous), yielding a transfer matrix at that interface. All these transfer and 
propagation matrices can then be combined to compute scattering/reflection 
matrices for an entire structure. 

The main difference here from classical RCWA computations is that the 
modes are computed at imaginary frequencies ^. As in Sec. |7.5.01 this actu- 
ally simplifies the problem. At an imaginary frequency cu = z^, the modes of 
a given cross-section e{i(,, x, y) and /i(i^, x, y) with z dependence e*^^ = e~^^ 
(7 = —i(3) satisfy the eigenequation (for isotropic materials) [134 135 



^^ See also the chapter by A. Lambrecht et al. in this volume for additional discussion of 
Casimir interactions among periodic structures. 
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, (7.26) 



where the xy subscript indicates a two-component vector with xy (trans- 
verse) components. The operators on both the left- and right-hand sides are 
real-symmetric, while the operator on the left-hand side is positive-definite, 
and as a result the eigenvalues 7 are purely real. This means that the prop- 
agation constants j3 are purely imaginary (all of the imaginary-frequency 
modes are evanescent in z), and the analogues of incoming/outgoing waves 
are those that are exponentially decaying towards/away from the surface. 
Moreover, the numerical problem of solving for these eigenmodes in a given 
cross-section reduces to a positive-definite generalized eigenvalue problem (a 
definite matrix pencil |69 ), to which the most desirable numerical solvers ap- 
ply (68l[69| (unlike the classical real-cj problem in which there are both prop- 
agating and evanescent modes because the problem is indefinite |134,135 ). 
For homogeneous cross-sections (as in the space between the two objects), 
the solutions are simply planewaves of the form ^'^^xx+ikyy-^fz+e^t ^ where for 
vacuum 7 = iy^lk^^^p + ^^/c^- 

For sufficiently simple cross-sections, especially in two-dimensional or ax- 
isymmetric geometries, it is possible to solve for the modes analytically and 
hence obtain the scattering matrices, and this is how the technique was first 
applied to the Casimir problem |18 . For more general geometries, one can 
solve for the modes numerically by a variety of techniques, such as by a 
transfer-matrix method in two dimensions |43 or by a planewave expansion 
(in the xy cross-section) in three dimensions |46 . Of course, one truncates to a 
finite number of modes via some cutoff I7I (which follows automatically from 
discretizing the cross-section in a finite grid, for example), and convergence is 
obtained in the limit as this cutoff increases. Given a basis of eigenmodes with 
some cutoff, the process of constructing the scattering/reffection matrices is 
thoroughly discussed elsewhere |43-47 , so we do not review it here. 

The strength of RCWA is that regions of uniform cross-section are han- 
dled with at most a 2d discretization of the cross-section, independent of the 
thickness of the region, so very thick or very thin layers can be solved effi- 
ciently. The main limitation of RCWA methods is that the transfer matrices 
(and the resulting reffection matrices Ri and R2) are dense N x N matrices, 
where N is the number of modes required for convergence. If N is large, as in 
complicated three-dimensional structures, the problem can quickly become 
impractical because of the 0{N'^) storage and 0{N^) computation require- 
ments. The most favorable case is that of periodic structures with relatively 
simple unit cells, in which case the problem can be reduced to that of com- 
puting the modes of each periodic unit cell (with Bloch-periodic boundary 
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conditions) as discussed below, and RCWA can then be quite practical even 
in three dimensions. Non-periodic structures, such as compact objects, can be 
handled by perfectly matched layer (PML) absorbing boundaries |44 , albeit 
at greater computational cost because of the increased cross-section size. 



7.7 Periodicity and other symmetries 

In this section, we briefly discuss the issue of periodicity and other symme- 
tries, which can be exploited to greatly reduce the computational effort in 
Casimir calculations just as for classical EM calculations. 

If a structure is periodic in the x direction with period yl, as in Fig. |7.5[ 
the problem simplifies considerably because one can reduce the computation 
to a single unit cell of thickness A. In particular, one imposes Bloch-periodic 
boundary conditions — the fields at x = yl equal the fields at x = mul- 
tiplied by a phase factor e*^=^^ — and computes the Casimir energy or force 
for each Bloch wavevector kx separately, then integrates the result over kx 
via j_^ij^{- • • )dkx- This can be derived in a variety of ways, for example by 
applying Bloch's theorem [g] to decompose the eigenmodes into Bloch-wave 
solutions for each kx^ or by expanding the delta functions of the fluctuation- 
dissipation approach in a Fourier series |12 . More generally, for any periodic 
unit cell, one can perform the Casimir energy/force computation for the unit 
cell with Bloch periodic boundaries and then integrate the Bloch wavevector 
k over the irreducible Brillouin zone (multiplied by the volume ratio of the 
Brillouin zone and the irreducible Brillouin zone). 

The specific case of continuous translational symmetry, say in the x di- 
rection, corresponds yl ^ and one must integrate over all kx (the Bril- 
louin zone is infinite). Certain additional simplifications apply in the case of 
a perfect-metal structure with continuous translational symmetry, in which 
case the fields decompose into two polarizations and the k integration can be 
performed implicitly [12]. 

Rotational symmetry can be handled similarly: the fields can be decom- 
posed into fields with e^^^ angular dependence, and the total force or energy 
is the sum over all integers m of the contributions for each m |27|. More gen- 
erally, the Casimir contributions can be decomposed into a sum of contribu- 
tions from irreducible representations of the symmetry group of the structure 
(e.g. all eigenmodes can be classified into these representations [136, 141| ); 



translational and rotational symmetries are merely special cases. As another 
example, in a structure with a mirror symmetry one could sum even- and 
odd-symmetry contributions (in fact, this is the underlying reason for the 
TE/TM polarization decomposition in two dimensions [6]). 
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7.8 Nonzero-temperature corrections 

In the preceding sections, we discussed only the computation of Casimir in- 
teractions at zero temperature T = 0+. However, the modification of any 
imaginary-frequency expression for a Casimir interaction from T = to T > 
is almost trivial: one simply performs a sum instead of an integral. If the 
T = interaction (energy, force, etc.) is expressed as an integral Jq C{^) d£^ 
of some contributions C(^) at each imaginary frequency ^, then the T > 
interaction is well known to be simply |81 : 

CX) 

Cm^^E'c(^n), (7.27) 

n=0 ^ ^ 

where /cb is Boltzmann's constant and ^ indicates a sum with weight ^ for 
the n = term. The frequencies ^^ = 2T:k^Tn/h are known as Matsubara 
frequencies, and the corresponding (imaginary) Matsubara wavelengths are 
\n = 27r/^n = At/^ where At = h/k^T. The conversion of the T = in- 
tegral into a summation can be derived in a variety of ways, most directly 
by considering thermodynamics in the Matsubara formalism |81 . Physically, 
this arises from the coXhihw / 2kT) Bose-Einstein distribution factor that ap- 



pears in the fluctuation-dissipation expressions (7.3) for nonzero tempera- 
tures. When the contour integration is performed over cj, the coth introduces 
poles at hw/2kT = inn that convert the integral into a sum via the residue 
theorem (with the n = residue having half weight because it lies on the 



real-cj axis) |77 . As explained in Sec. 7.5.1.2, some care must be applied in 



evaluating the n = term because of the well known singularity of Maxwell's 
equations at cj = (where the E and H fields decouple), and one may need to 
take the limit ^ ^ 0+ (although there is some controversy in the unphysical 
case of perfect metals |92 ) 



Mathematically, the sum of eq. (7.27) is exactly the same as a trapezoidal- 



rule approximation for t he T = integral, with equally spaced abscissas 
A^ = 27t/Xt |53,91[[137). Thanks to the 0(/A^^) convergence of the trape- 
zoidal rule |53 , this means that the T > result is quite close to the T = 
result unless C(^) varies rapidly on the scale of 27r/AT. In particular, sup- 
pose that C(^) varies on a scale 27r/a , corresponding to some lengthscale 
a in the problem (typically from a surface-surface separation). In that case, 
assuming C(^) has nonzero slopa^ at ^ = 0+ (typical for interactions be- 
tween realistic metal surfaces), then the nonzero-T correction should be of 
order 0(a^/A|.). At room temperature (T = 300 K), At ^ 7.6 /im, and the 
temperature corrections to Casimir interactions are typically negligible for 
submicron separations 1 81, 1381. On the other hand, it is possible that careful 
material and geometry choices may lead to larger temperature effects |137|. 



16 j£ C(^) has zero slope at ^ = 0+, then the trapezoidal rule differs from the integral by 
0(A^'^) or less, depending upon which derivative is nonzero at ^ = 0+ j53|. 
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There is also the possibihty of interesting effects in nonequihbrium situations 
(objects at different temperatures) |97{|139| , but such situations are beyond 
the scope of this review. 



7.9 Concluding remarks 

The area of numerical Casimir computations remains rich with opportuni- 
ties. Relatively few geometry and material combinations have as yet been 
explored, and thus many newly answerable questions remain regarding the 
ways in which Casimir phenomena can be modified by exploiting the degrees 
of freedom available in modern nanofabrication. In the regime of computa- 
tional techniques, while several effective methods have already been proposed 
and demonstrated, the parallels with computational electromagnetism lead 
us to anticipate ongoing improvements and developments for some time to 
come. The same parallels also caution against any absolute "rankings" of the 
different approaches, as different numerical techniques have always exhibited 
unique strengths and weaknesses in both theory and practice. And because 
computer time is typically much less expensive than programmer time, there 
is something to be said for methods that may be theoretically suboptimal 
but are easy to implement (or are available off-the-shelf) for very general 
geometries and materials. Nor is the value of analytical and semi- analytical 
techniques diminished, but rather these approaches are freed from the tedium 
of hand computation to focus on more fundamental questions. 
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